Lab 1 2021
Mark Hebblewhite
1/16/2021
1 Lab Objectives
• Introduce students to R through and R studio, a Graphic User Interface (GUI) • REVIEW basic principle of data management, data exploration and visualization, and basic statistical tests including t-tests, linear models, and regression. • Introduce the concept of scripting to conduct statistical analyses. • Introduce how to handle spatial data (Rasters, Vectors) using the R packages sp, raster.
2 1.0 Why R? R & R Studio
R is a free, open-source statistics software program. R is a command-line driven high-level software language, meaning a user enters code into the R console to make calculations, conduct analyses, etc. Unfortunately, this style of software can be daunting for beginners, who are not familiar with the commands or syntax of R – and because R is free the help files can sometimes be less than helpful, but they are improving. R Commander smoothes much of this confusion and aggravation by providing a point-and-click window environment for new users to begin performing basic analyses. In addition to these user-friendly windows, you can also enter R code directly into the script window in R commander and click the submit button to run it. The script window serves as the R console. Note, however, this friendly environment may also sacrifice some of the flexibility of using the full R language and may become frustrating as you learn the R language. In summary, there are 3 main reasons why the ecological world is moving to R; 1) It is free, and has a growing and loyal on-line community of software developers continually improving and adding to its functionality – chances are, if you have to do a customized analysis, someone has already built a package to download and use. 2) Its becoming the statistical software of choice of ecologists and statisticians worldwide, and will carry your career for at least the next decade. 3) Most importantly, though, it is a fully customizable software platform (not a programming language) for conducting any analysis you need in your graduate research.
R Resources • Website - http://www.r-project.org/index.html • The R-journal - http://journal.r-project.org/ • Internet searching - http://www.rseek.org/
R studio https://www.rstudio.com Another of Hadley Wickham’s creations (along with ggplot2, Rmarkdown, etc), R studio is a RStudio is an integrated development environment (IDE) for R.
3 Preliminaries
3.1 Find out about R
R.Version()## $platform
## [1] "x86_64-w64-mingw32"
##
## $arch
## [1] "x86_64"
##
## $os
## [1] "mingw32"
##
## $crt
## [1] "ucrt"
##
## $system
## [1] "x86_64, mingw32"
##
## $status
## [1] ""
##
## $major
## [1] "4"
##
## $minor
## [1] "2.2"
##
## $year
## [1] "2022"
##
## $month
## [1] "10"
##
## $day
## [1] "31"
##
## $`svn rev`
## [1] "83211"
##
## $language
## [1] "R"
##
## $version.string
## [1] "R version 4.2.2 (2022-10-31 ucrt)"
##
## $nickname
## [1] "Innocent and Trusting"citation()##
## To cite R in publications use:
##
## R Core Team (2022). R: A language and environment for statistical
## computing. R Foundation for Statistical Computing, Vienna, Austria.
## URL https://www.R-project.org/.
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {R: A Language and Environment for Statistical Computing},
## author = {{R Core Team}},
## organization = {R Foundation for Statistical Computing},
## address = {Vienna, Austria},
## year = {2022},
## url = {https://www.R-project.org/},
## }
##
## We have invested a lot of time and effort in creating R, please cite it
## when using it for data analysis. See also 'citation("pkgname")' for
## citing R packages.4 Objective 1.0 Cougar Track Count Analysis
First we will refresh our data exploration, summary statistic skills, and generalized linear modeling skills to review some basic foundations of statistics and R together.
We will be using data on track counts of mountain lions/cougars in and around Banff National Park Alberta in the data file called cougar.csv
The file cougar.csv contains data on cougar track counts from the Bow Valley of Banff National Park. Over 5 years, multiple cougars were snow tracked through various wildlife corridors surrounding the townsite of Banff (Duke et al., 2001), and the number of counts of tracks of cougars within 30m2 pixels was counted in one wildlife corridor. This resulted in some cells having multiple counts of cougar tracks and others having no counts. Park Managers wanted to know what habitat factors cougars needed for corridor use, and the dataset includes information about the count of cougar tracks and associated spatial covariates including slope, distance to cover, distance to roads and distance to trails.
Duke, D. L., M. Hebblewhite, P. C. Paquet, C. Callaghan, and M. Percy.
2001. Restoration of a large carnivore corridor in Banff National Park.
Pages 261-275 in D. S. Maehr, R. F. Noss, and J. L. Larkin, editors.
Large mammal restoration: ecological and sociological challenges in the
21st century. Island Press, Washington.
Import the data csv
cougar <- read.csv("Data/cougar.csv", header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE)
head(cougar)## Use UseNonUse Slope CoverDist Roads AllTrails
## 1 0 0 9.462322 0 787 228
## 2 0 0 39.856209 0 690 247
## 3 0 0 0.000000 7 150 153
## 4 0 0 11.386165 0 30 671
## 5 0 0 10.555381 0 457 42
## 6 0 0 7.653049 0 335 192str(cougar)## 'data.frame': 726 obs. of 6 variables:
## $ Use : int 0 0 0 0 0 0 0 0 0 0 ...
## $ UseNonUse: int 0 0 0 0 0 0 0 0 0 0 ...
## $ Slope : num 9.46 39.86 0 11.39 10.56 ...
## $ CoverDist: int 0 0 7 0 0 0 0 0 16 0 ...
## $ Roads : int 787 690 150 30 457 335 741 711 630 283 ...
## $ AllTrails: int 228 247 153 671 42 192 60 42 300 134 ...table(cougar$Use, cougar$UseNonUse)##
## 0 1
## 0 363 0
## 1 0 31
## 2 0 34
## 3 0 28
## 4 0 32
## 5 0 28
## 6 0 121
## 7 0 16
## 8 0 14
## 9 0 10
## 10 0 12
## 11 0 8
## 12 0 18
## 13 0 4
## 14 0 5
## 15 0 1
## 19 0 1# Cougar data management
## Let us convert the USENONUSE column to a factor for future analysis, and store it as as a factorUSE
cougar$factorUSE <- as.factor(cougar$UseNonUse)4.1 Cougar summary statistics
Lets calculate some basic summary statistics for the data set…
summary(cougar)## Use UseNonUse Slope CoverDist
## Min. : 0.000 Min. :0.0 Min. : 0.000 Min. : 0.000
## 1st Qu.: 0.000 1st Qu.:0.0 1st Qu.: 0.000 1st Qu.: 0.000
## Median : 0.500 Median :0.5 Median : 2.134 Median : 0.000
## Mean : 2.853 Mean :0.5 Mean : 8.111 Mean : 2.959
## 3rd Qu.: 6.000 3rd Qu.:1.0 3rd Qu.:12.773 3rd Qu.: 0.000
## Max. :19.000 Max. :1.0 Max. :47.730 Max. :53.000
## Roads AllTrails factorUSE
## Min. : 0.0 Min. : 0.0 0:363
## 1st Qu.: 260.5 1st Qu.: 67.0 1:363
## Median : 450.0 Median :175.0
## Mean : 465.5 Mean :224.8
## 3rd Qu.: 641.0 3rd Qu.:323.0
## Max. :1400.0 Max. :899.0# ... by used and unused locations
install.packages("tidyverse")## package 'tidyverse' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\Administrator.KJLWS11\AppData\Local\Temp\Rtmp84rxPr\downloaded_packageslibrary(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0 ✔ purrr 1.0.0
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.5.0
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()cougar_df <- as_tibble(cougar)
byUse <- group_by(cougar_df, UseNonUse)
summarise(byUse, slope = mean(Slope))## # A tibble: 2 × 2
## UseNonUse slope
## <int> <dbl>
## 1 0 10.3
## 2 1 5.93summarise(byUse, DistTrails = mean(AllTrails))## # A tibble: 2 × 2
## UseNonUse DistTrails
## <int> <dbl>
## 1 0 244.
## 2 1 206.summarise(byUse, DistCover = mean(CoverDist))## # A tibble: 2 × 2
## UseNonUse DistCover
## <int> <dbl>
## 1 0 3.84
## 2 1 2.07summarise(byUse, DistRoads = mean(Roads))## # A tibble: 2 × 2
## UseNonUse DistRoads
## <int> <dbl>
## 1 0 474.
## 2 1 457.4.2 Cougar graphing
Lets visualize this differences with overlaid histograms for each variable. The multhist function from the plotrix package does a nice job of doing this.
install.packages("plotrix")## package 'plotrix' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\Administrator.KJLWS11\AppData\Local\Temp\Rtmp84rxPr\downloaded_packageslibrary(plotrix)
# Use a 2x2 plotting matrix to see all of the histograms at once
par(mfrow = c(2,2))
multhist(list(cougar$AllTrails[cougar$factorUSE==1],cougar$AllTrails[cougar$factorUSE==0]), freq = TRUE, main = "Trails")
# I chose to put a legend in the upper right hand graph.
# That's what the additional arguments in the line below specify.
multhist(list(cougar$CoverDist[cougar$factorUSE==1],cougar$CoverDist[cougar$factorUSE==0]), freq = TRUE, main = "Cover Distance", legend.text = c("Used", "Unused"), args.legend = list(bty = "n"))
multhist(list(cougar$Roads[cougar$factorUSE==1],cougar$Roads[cougar$factorUSE==0]), freq = TRUE, main = "Roads")
multhist(list(cougar$Slope[cougar$factorUSE==1],cougar$Slope[cougar$factorUSE==0]), freq = TRUE, main = "Slope")
4.2.1 Cougar Boxplots
Box plots also do a good job of visualizing these differences. These are plotted in a similar fashion as above.
par(mfrow= c(2,2))
boxplot(AllTrails~factorUSE, ylab="Distance (m)", xlab="Used",main = "Trails", data=cougar)
boxplot(CoverDist~factorUSE, ylab="Distance (m)", xlab="Used", main = "Cover", data=cougar)
boxplot(Roads~factorUSE, ylab="Distance (m)", xlab="Used",main = "Roads", data=cougar)
boxplot(Slope~factorUSE, ylab="Slope", xlab="Used", main = "Slope", data=cougar)
4.3 Cougar statistical tests
Test for significant differences in each of the covariates
t.test(AllTrails~factorUSE, alternative='two.sided', conf.level=.95,
var.equal=FALSE, data=cougar)##
## Welch Two Sample t-test
##
## data: AllTrails by factorUSE
## t = 2.6896, df = 715.58, p-value = 0.00732
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## 10.13185 64.90396
## sample estimates:
## mean in group 0 mean in group 1
## 243.5399 206.0220t.test(CoverDist~factorUSE, alternative='two.sided', conf.level=.95,
var.equal=FALSE, data=cougar)##
## Welch Two Sample t-test
##
## data: CoverDist by factorUSE
## t = 3.3686, df = 604.85, p-value = 0.0008037
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## 0.7374901 2.7997000
## sample estimates:
## mean in group 0 mean in group 1
## 3.842975 2.074380t.test(Roads~factorUSE, alternative='two.sided', conf.level=.95,
var.equal=FALSE, data=cougar)##
## Welch Two Sample t-test
##
## data: Roads by factorUSE
## t = 0.86435, df = 679.54, p-value = 0.3877
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -22.24071 57.22143
## sample estimates:
## mean in group 0 mean in group 1
## 474.2039 456.7135t.test(Slope~factorUSE, alternative='two.sided', conf.level=.95,
var.equal=FALSE, data=cougar)##
## Welch Two Sample t-test
##
## data: Slope by factorUSE
## t = 5.4158, df = 627.87, p-value = 8.7e-08
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## 2.774890 5.931965
## sample estimates:
## mean in group 0 mean in group 1
## 10.287312 5.9338844.4 Cougar statistical analyses - linear models
First, letes create a data set of only the used locations
USEonly <- subset(cougar, subset=UseNonUse == 1)
# Construct linear models of track counts as a function of 4 covariates
# Distance to trails model
trails <- glm(Use ~ AllTrails, family=gaussian(identity), data=USEonly)
summary(trails)##
## Call:
## glm(formula = Use ~ AllTrails, family = gaussian(identity), data = USEonly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.120 -2.106 -0.033 1.054 13.009
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.1910009 0.2541896 24.356 <2e-16 ***
## AllTrails -0.0023578 0.0009354 -2.521 0.0121 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 9.972573)
##
## Null deviance: 3663.5 on 362 degrees of freedom
## Residual deviance: 3600.1 on 361 degrees of freedom
## AIC: 1869
##
## Number of Fisher Scoring iterations: 2# Slope model
slope <- glm(Use ~ Slope, family=gaussian(identity), data=USEonly)
summary(slope)##
## Call:
## glm(formula = Use ~ Slope, family = gaussian(identity), data = USEonly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.7694 -2.6853 0.2811 0.3147 13.3114
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.685274 0.204401 27.81 <2e-16 ***
## Slope 0.003364 0.019816 0.17 0.865
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 10.14728)
##
## Null deviance: 3663.5 on 362 degrees of freedom
## Residual deviance: 3663.2 on 361 degrees of freedom
## AIC: 1875.3
##
## Number of Fisher Scoring iterations: 2# Distance to cover model
cover <- glm(Use ~ CoverDist, family=gaussian(identity), data=USEonly)
summary(cover)##
## Call:
## glm(formula = Use ~ CoverDist, family = gaussian(identity), data = USEonly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8123 -1.9863 0.1877 1.0137 13.1877
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.81233 0.17904 32.464 <2e-16 ***
## CoverDist -0.05163 0.03162 -1.633 0.103
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 10.07372)
##
## Null deviance: 3663.5 on 362 degrees of freedom
## Residual deviance: 3636.6 on 361 degrees of freedom
## AIC: 1872.6
##
## Number of Fisher Scoring iterations: 2# Distance to roads model
roads <- glm(Use ~ Roads, family=gaussian(identity), data=USEonly)
summary(roads)##
## Call:
## glm(formula = Use ~ Roads, family = gaussian(identity), data = USEonly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8972 -2.5408 0.2199 0.5703 13.2917
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.9166147 0.3654078 16.19 <2e-16 ***
## Roads -0.0004628 0.0007115 -0.65 0.516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 10.13621)
##
## Null deviance: 3663.5 on 362 degrees of freedom
## Residual deviance: 3659.2 on 361 degrees of freedom
## AIC: 1874.9
##
## Number of Fisher Scoring iterations: 2But these tests require the assumption of normality for the response variable, use. Is that reasonable? Lets test track counts for normality
# Visualize with a histogram
par(mfrow= c(1,1))
hist(USEonly$Use, scale="frequency", breaks="Sturges", col="darkgray")## Warning in plot.window(xlim, ylim, "", ...): "scale" is not a graphical
## parameter## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...): "scale"
## is not a graphical parameter## Warning in axis(1, ...): "scale" is not a graphical parameter## Warning in axis(2, at = yt, ...): "scale" is not a graphical parameter
## not really that normal
shapiro.test(USEonly$Use)##
## Shapiro-Wilk normality test
##
## data: USEonly$Use
## W = 0.92147, p-value = 7.473e-13# Visualize with a histogramclearly, the answer is NO. This means we need to consider OTHER forms of generalized linear models, or transformation of the response variable, here, USE.
Lets re-fit linear models with a ln transform, first make the ln transform, and then refit the linear models with the transformed use
USEonly$lnUSE <- with(USEonly, log(Use))
hist(USEonly$lnUSE) ## a bit more normal
### Now re-fit the models
# Distance to trails model
ln.trails <- glm(lnUSE ~ AllTrails, family=gaussian(identity), data=USEonly)
summary(ln.trails)##
## Call:
## glm(formula = lnUSE ~ AllTrails, family = gaussian(identity),
## data = USEonly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6225 -0.2778 0.1812 0.3611 1.3439
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.6345317 0.0541846 30.166 <2e-16 ***
## AllTrails -0.0003999 0.0001994 -2.006 0.0456 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.4531516)
##
## Null deviance: 165.41 on 362 degrees of freedom
## Residual deviance: 163.59 on 361 degrees of freedom
## AIC: 746.82
##
## Number of Fisher Scoring iterations: 2# Slope model
ln.slope <- glm(lnUSE ~ Slope, family=gaussian(identity), data=USEonly)
summary(ln.slope)##
## Call:
## glm(formula = lnUSE ~ Slope, family = gaussian(identity), data = USEonly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5607 -0.4340 0.2311 0.2707 1.3852
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.56068 0.04343 35.938 <2e-16 ***
## Slope -0.00144 0.00421 -0.342 0.732
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.4580528)
##
## Null deviance: 165.41 on 362 degrees of freedom
## Residual deviance: 165.36 on 361 degrees of freedom
## AIC: 750.72
##
## Number of Fisher Scoring iterations: 2# Distance to cover model
ln.cover <- glm(lnUSE ~ CoverDist, family=gaussian(identity), data=USEonly)
summary(ln.cover)##
## Call:
## glm(formula = lnUSE ~ CoverDist, family = gaussian(identity),
## data = USEonly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5792 -0.2721 0.2126 0.3667 1.4028
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.57917 0.03799 41.573 <2e-16 ***
## CoverDist -0.01303 0.00671 -1.942 0.0529 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.4534625)
##
## Null deviance: 165.41 on 362 degrees of freedom
## Residual deviance: 163.70 on 361 degrees of freedom
## AIC: 747.07
##
## Number of Fisher Scoring iterations: 2# Distance to roads model
ln.roads <- glm(lnUSE ~ Roads, family=gaussian(identity), data=USEonly)
summary(ln.roads)##
## Call:
## glm(formula = lnUSE ~ Roads, family = gaussian(identity), data = USEonly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5737 -0.4394 0.2330 0.2563 1.3926
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.535e+00 7.768e-02 19.76 <2e-16 ***
## Roads 3.781e-05 1.513e-04 0.25 0.803
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.4581219)
##
## Null deviance: 165.41 on 362 degrees of freedom
## Residual deviance: 165.38 on 361 degrees of freedom
## AIC: 750.78
##
## Number of Fisher Scoring iterations: 2Next, lets switch back to the complete data set, and then re-do the Shapiro-Wilk Test and make another histogram for all of the data.
shapiro.test(cougar$Use)##
## Shapiro-Wilk normality test
##
## data: cougar$Use
## W = 0.78045, p-value < 2.2e-16hist(cougar$Use, scale="frequency", breaks="Sturges", col="darkgray")## Warning in plot.window(xlim, ylim, "", ...): "scale" is not a graphical
## parameter## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...): "scale"
## is not a graphical parameter## Warning in axis(1, ...): "scale" is not a graphical parameter## Warning in axis(2, at = yt, ...): "scale" is not a graphical parameter
Clearly these data are still not normal, and so it might be useful to
think of categorizing these data as USED or NOT USED, in a binary
fashion, and using logistic regression of the use = 1 and not use = 0.
Thus, next we will use logistic regression to model the log-odds of USE
for each of the 4 covariates
# Trails model
logitTrails <- glm(UseNonUse ~ AllTrails, family=binomial(logit), data=cougar)
summary(logitTrails)##
## Call:
## glm(formula = UseNonUse ~ AllTrails, family = binomial(logit),
## data = cougar)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2805 -1.1920 0.0936 1.1450 1.4275
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2392204 0.1163257 2.056 0.03974 *
## AllTrails -0.0010675 0.0004007 -2.664 0.00772 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1006.45 on 725 degrees of freedom
## Residual deviance: 999.22 on 724 degrees of freedom
## AIC: 1003.2
##
## Number of Fisher Scoring iterations: 4# Slope model
logitSlope <- glm(UseNonUse ~ Slope, family=binomial(logit), data=cougar)
summary(logitSlope)##
## Call:
## glm(formula = UseNonUse ~ Slope, family = binomial(logit), data = cougar)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3078 -1.2388 0.2123 1.0679 1.8429
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.301367 0.093849 3.211 0.00132 **
## Slope -0.038241 0.007427 -5.149 2.62e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1006.45 on 725 degrees of freedom
## Residual deviance: 977.23 on 724 degrees of freedom
## AIC: 981.23
##
## Number of Fisher Scoring iterations: 4# Cover model
logitCover <- glm(UseNonUse~ CoverDist, family=binomial(logit), data=cougar)
summary(logitCover)##
## Call:
## glm(formula = UseNonUse ~ CoverDist, family = binomial(logit),
## data = cougar)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2233 -1.2233 0.3021 1.1322 1.8036
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.10739 0.08089 1.328 0.18428
## CoverDist -0.03788 0.01175 -3.224 0.00126 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1006.45 on 725 degrees of freedom
## Residual deviance: 994.78 on 724 degrees of freedom
## AIC: 998.78
##
## Number of Fisher Scoring iterations: 4# Roads model
logitRoads <- glm(UseNonUse ~ Roads, family=binomial(logit), data=cougar)
summary(logitRoads)##
## Call:
## glm(formula = UseNonUse ~ Roads, family = binomial(logit), data = cougar)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.22439 -1.18003 0.02292 1.17627 1.24756
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1098329 0.1471202 0.747 0.455
## Roads -0.0002360 0.0002729 -0.865 0.387
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1006.4 on 725 degrees of freedom
## Residual deviance: 1005.7 on 724 degrees of freedom
## AIC: 1009.7
##
## Number of Fisher Scoring iterations: 34.5 Logistic Regression Visualization - graphing
install.packages("ggplot2")## Warning: package 'ggplot2' is in use and will not be installedlibrary(ggplot2)
ggplot(cougar, aes(x=Slope, y=UseNonUse)) + geom_rug() + stat_smooth(method="glm", method.args=list(family="binomial"))## `geom_smooth()` using formula = 'y ~ x'
ggplot(cougar, aes(x=CoverDist, y=UseNonUse)) + geom_rug() + stat_smooth(method="glm", method.args=list(family="binomial"))## `geom_smooth()` using formula = 'y ~ x'
So, to review, we have reviewed some basic data visualization steps (boxplots, histograms, etc), basic summary statistics, basic statistical analyses (t-tests, tests for normality, generalized linear models), and considered a simple dataset with different ways of thinking about animal use of space and 4 spatial covariates.
But, how did we build our dataset of things like distance to cover, distance to trails, etc? Next, we will move to understanding and exploring spatial data in R.
5 Objective 2 - exploring SPATIAL data in R.
This semester we will use this nifty little wrapper package Dan Eacker adapted from here: https://gist.github.com/stevenworthington/3178163 ## Define function to install and load required packages
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if (length(new.pkg))
install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
#load or install these packages:
packages <- c("ks", "lattice", "adehabitatHR", "maptools", "foreign", "rgdal", "sp", "raster","plot3D","rasterVis", "colorRamps","rgeos")
#run function to install packages - e.g., library command
ipak(packages)## Loading required package: ks## Loading required package: lattice## Loading required package: adehabitatHR## Loading required package: sp## Loading required package: deldir## deldir 1.0-6 Nickname: "Mendacious Cosmonaut"##
## The syntax of deldir() has had an important change.
## The arguments have been re-ordered (the first three
## are now "x, y, z") and some arguments have been
## eliminated. The handling of the z ("tags")
## argument has been improved.
##
## The "dummy points" facility has been removed.
## This facility was a historical artefact, was really
## of no use to anyone, and had hung around much too
## long. Since there are no longer any "dummy points",
## the structure of the value returned by deldir() has
## changed slightly. The arguments of plot.deldir()
## have been adjusted accordingly; e.g. the character
## string "wpoints" ("which points") has been
## replaced by the logical scalar "showpoints".
## The user should consult the help files.## Loading required package: ade4## Loading required package: adehabitatMA## Registered S3 methods overwritten by 'adehabitatMA':
## method from
## print.SpatialPixelsDataFrame sp
## print.SpatialPixels sp## Loading required package: adehabitatLT## Loading required package: CircStats## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: boot##
## Attaching package: 'boot'## The following object is masked from 'package:lattice':
##
## melanoma##
## Attaching package: 'adehabitatLT'## The following object is masked from 'package:dplyr':
##
## id## Loading required package: maptools## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired during 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.## Loading required package: foreign## Loading required package: rgdal## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-3, (SVN revision 1187)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Program Files/R/R-4.2.2/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Program Files/R/R-4.2.2/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.## Loading required package: raster##
## Attaching package: 'raster'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:dplyr':
##
## select## Loading required package: plot3D## Loading required package: rasterVis## Loading required package: colorRamps## Loading required package: rgeos## rgeos version: 0.6-1, (SVN revision 692)
## GEOS runtime version: 3.9.3-CAPI-1.14.3
## Please note that rgeos will be retired during 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
## GEOS using OverlayNG
## Linking to sp version: 1.5-1
## Polygon checking: TRUE## ks lattice adehabitatHR maptools foreign rgdal
## TRUE TRUE TRUE TRUE TRUE TRUE
## sp raster plot3D rasterVis colorRamps rgeos
## TRUE TRUE TRUE TRUE TRUE TRUE5.1 Working Directly with Shapefiles in R
First set working directory, then read in shapefiles
# reading in shapefiles (raster package)
elc_habitat<-shapefile("Data/elc_habitat.shp")
humanaccess<-shapefile("Data/humanacess.shp")
mcp2<-shapefile("Data/mcp2.shp")
wolfyht<-shapefile("Data/wolfyht.shp")
# make a very basic plot of shapefile after resetting graphical parameters
par(mfrow= c(1,1))
plot(elc_habitat)
plot(wolfyht)
# look at the class of the shapefile
class(elc_habitat)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"# look at structure of shapefile
#str(elc_habitat)
# look at first 20 rows of data for shapefile
head(elc_habitat@data, n=20)## ID ECO_NUM CODE ECOREGION ECOSECTION SLOPE COMPLEX MODIFIER_1
## 1 1 154 WF3 Upper Subalpine Wildflower 9 <NA> B
## 2 2 101 PL1 Upper Subalpine Peyto Lake 8 <NA> <NA>
## 3 3 152 WF1 Upper Subalpine Wildflower 9 <NA> A
## 4 4 125 SB2 Lower Subalpine Sawback 8 <NA> A
## 5 5 153 WF2 Upper Subalpine Wildflower 8 <NA> <NA>
## 6 6 110 PP6 Lower Subalpine Pipestone 7 <NA> B
## 7 7 127 SB4 Lower Subalpine Sawback 9 <NA> B
## 8 8 1 AL1 Lower Subalpine Altrude 3 <NA> B
## 9 9 152 WF1 Upper Subalpine Wildflower 8 <NA> A
## 10 10 70 HC4 Lower Subalpine Hector 3 c <NA>
## 11 11 108 PP3 Lower Subalpine Pipestone 7 <NA> <NA>
## 12 12 154 WF3 Upper Subalpine Wildflower 9 <NA> A
## 13 13 180 ZZ rock_water Water 0 <NA> <NA>
## 14 14 152 WF1 Upper Subalpine Wildflower 9 <NA> X
## 15 15 112 PR1 Lower Subalpine Panarama 7 c B
## 16 16 179 T Alpine Rock 0 <NA> <NA>
## 17 17 9 BK1 Lower Subalpine Baker_cr 6 c <NA>
## 18 18 153 WF2 Upper Subalpine Wildflower 8 <NA> <NA>
## 19 19 181 R+CR rock_water Rock 0 <NA> <NA>
## 20 20 180 ZZ rock_water Water 0 <NA> <NA>
## MODIFIER_2 DESCRIPTIO DEFINITION DESCRIPT0 EREG_NUM ESEC_NUM
## 1 <NA> BURNED > 70% <NA> 3 55
## 2 <NA> <NA> 45 - 70% <NA> 3 39
## 3 B SNOW AVALANCHED > 70% BURNED 3 55
## 4 <NA> SNOW AVALANCHED 45 - 70% <NA> 2 45
## 5 <NA> <NA> 45 - 70% <NA> 3 55
## 6 <NA> BURNED 30 - 45% <NA> 2 40
## 7 <NA> BURNED > 70% <NA> 2 45
## 8 <NA> BURNED 0 - 5% <NA> 2 1
## 9 <NA> SNOW AVALANCHED 45 - 70% <NA> 3 55
## 10 <NA> <NA> 0 - 5% <NA> 2 24
## 11 <NA> <NA> 30 - 45% <NA> 2 40
## 12 <NA> SNOW AVALANCHED > 70% <NA> 3 55
## 13 <NA> <NA> <NA> <NA> 7 59
## 14 <NA> LITHIC > 70% <NA> 3 55
## 15 <NA> BURNED 30 - 45% <NA> 2 41
## 16 <NA> <NA> <NA> <NA> 6 61
## 17 <NA> <NA> 15 - 30% <NA> 2 4
## 18 <NA> <NA> 45 - 70% <NA> 3 55
## 19 <NA> <NA> <NA> <NA> 7 61
## 20 <NA> <NA> <NA> <NA> 7 59
## DOMVEG DVEG_NUM VEGTYPE VTYP_NUM DOMSOIL
## 1 Open_Mixed_conifer 16 O4 53 Brunisol
## 2 Englemann_spruce/Subalpine_fir 3 C15 14 Brunisol
## 3 Englemann_spruce/Subalpine_fir 3 C15 14 Brunisol
## 4 Englemann_spruce/Subalpine_fir 3 C13 12 Brunisol
## 5 Salix_shrubland 22 S2 69 Brunisol
## 6 Open_Englemann_spruce/Subalpine_fir 4 O18 65 Regosolic
## 7 Lodgepole_pine 1 O4 53 Brunisol
## 8 Lodgepole_pine 1 C19 18 Non-soil
## 9 Englemann_spruce/Subalpine_fir 3 C15 14 Brunisol
## 10 Dwarf_birch_shrubland 5 S1 68 Gleysolic
## 11 Englemann_spruce/Subalpine_fir 3 C13 12 Regosolic
## 12 Open_Mixed_conifer 16 O4 53 Brunisol
## 13 Non_vegetated 29 NIL 111 Non-soil
## 14 Englemann_spruce/Subalpine_fir 3 C15 14 Brunisol
## 15 Lodgepole_pine 1 C20 19 Brunisol
## 16 Non_vegetated 29 NIL 111 Non-soil
## 17 Lodgepole_pine 1 C18 17 Brunisol
## 18 Salix_shrubland 22 S2 69 Brunisol
## 19 Non_vegetated 29 NIL 111 Non-soil
## 20 Non_vegetated 29 NIL 111 Non-soil
## DSOI_NUM CLS DCLS_NUM LANDFORM LFOR_NUM DEERHAB_W
## 1 1 <NA> 0 Colluvial 2 0
## 2 1 Moderately 4 Till 12 3
## 3 1 <NA> 0 Colluvial 2 0
## 4 1 <NA> 0 Colluvial 2 0
## 5 1 <NA> 0 Colluvial 2 0
## 6 6 Well 3 Fluvial 4 3
## 7 1 <NA> 0 Colluvial 2 0
## 8 1 Well 3 Fluvial 4 3
## 9 1 <NA> 0 Colluvial 2 0
## 10 5 Poorly 6 Fluvial 4 3
## 11 6 Well 3 Fluvial 4 3
## 12 1 <NA> 0 Colluvial 2 0
## 13 8 <NA> 0 Water 0 0
## 14 1 <NA> 0 Colluvial 2 0
## 15 1 Well 3 Till 12 4
## 16 8 <NA> 0 Talus 2 0
## 17 1 Well 3 morainal 8 3
## 18 1 <NA> 0 Colluvial 2 0
## 19 8 <NA> 0 Rock and Colluvial Rubble 13 0
## 20 8 <NA> 0 Water 0 0
## DEERHAB_S E00_CODE CODE0 VEGROUP CARI_W CARI_S DEER_W DEER_S MOOSE_W
## 1 0 011248WF3B/9 WF3B <NA> 1 1 3 1 3
## 2 2 011114PL1/8 PL1 <NA> 5 3 4 3 4
## 3 0 011223WF1AB/9 WF1AB <NA> 1 1 1 1 4
## 4 0 011198SB2A/8 SB2A <NA> 1 1 3 1 4
## 5 0 011230WF2/8 WF2 <NA> 3 5 4 3 3
## 6 2 011164PP6B/7 PP6B <NA> 4 4 4 3 5
## 7 0 011209SB4B/9 SB4B <NA> 1 1 4 1 3
## 8 4 011008AL1B/3 AL1 <NA> 1 1 4 5 4
## 9 0 011221WF1A/8 WF1A <NA> 1 1 1 1 4
## 10 3 011085HC4/3c HC4 <NA> 5 3 4 4 5
## 11 2 011160PP3/7 PP3 <NA> 3 1 4 3 5
## 12 0 011242WF3A/9 WF3A <NA> 1 1 3 1 3
## 13 0 011260ZZ ZZ <NA> 7 7 7 7 7
## 14 0 011226WF1X/9 WF1X <NA> 1 1 1 1 4
## 15 2 011172PR1B/7c PR1B <NA> 1 1 5 3 4
## 16 0 011214T T <NA> 1 1 1 1 1
## 17 4 011017BK1/6c BK1 <NA> 1 1 4 5 5
## 18 0 011230WF2/8 WF2 <NA> 3 5 4 3 3
## 19 0 011188R+CR R+CR <NA> 1 1 1 1 1
## 20 0 011260ZZ ZZ <NA> 7 7 7 7 7
## MOOSE_S ELK_W ELK_S GOAT_W GOAT_S SHEEP_W SHEEP_S WOLF_W WOLF_S COYOTE_W
## 1 4 3 3 6 5 5 5 3 3 3
## 2 4 4 4 3 3 4 4 4 4 4
## 3 4 1 4 4 1 3 3 3 4 3
## 4 4 4 4 4 1 3 4 4 3 4
## 5 4 3 3 6 5 5 5 4 4 4
## 6 4 5 5 1 3 3 3 5 5 5
## 7 1 4 3 5 4 5 4 4 3 4
## 8 3 4 5 3 1 1 1 4 5 4
## 9 4 1 4 4 1 3 3 3 4 3
## 10 5 5 5 3 3 3 3 5 5 5
## 11 4 4 4 3 3 3 3 4 4 4
## 12 4 3 3 6 5 5 5 3 3 3
## 13 7 7 7 7 7 7 7 7 7 7
## 14 4 1 4 4 1 3 3 3 4 3
## 15 3 3 3 3 1 3 1 5 3 4
## 16 1 1 1 3 3 1 1 1 1 1
## 17 5 4 5 1 1 3 3 4 5 5
## 18 4 3 3 6 5 5 5 4 4 4
## 19 1 1 1 4 3 1 1 1 1 1
## 20 7 7 7 7 7 7 7 7 7 7
## COYOTE_S MARTEN FISHER WEASEL MINK WOLVERINE COUGAR_W COUGAR_S LYNX PIKA
## 1 3 1 1 1 1 1 3 3 3 4
## 2 4 5 1 5 1 4 3 4 3 5
## 3 4 5 1 1 1 1 1 3 1 4
## 4 3 3 1 1 1 1 4 4 4 4
## 5 4 3 1 3 1 5 3 4 1 5
## 6 5 4 1 5 1 4 5 4 5 1
## 7 3 1 1 1 1 1 5 3 4 5
## 8 5 4 1 4 1 5 4 5 5 1
## 9 4 5 1 1 1 1 1 3 1 4
## 10 6 4 2 5 2 3 5 5 4 3
## 11 4 4 1 4 1 4 4 4 5 3
## 12 3 1 1 1 1 1 3 3 3 4
## 13 7 7 7 7 7 7 7 7 7 7
## 14 4 5 1 1 1 1 1 3 1 4
## 15 4 4 2 3 2 5 4 3 4 4
## 16 1 1 1 4 1 1 1 1 1 1
## 17 6 5 1 3 2 5 4 5 4 4
## 18 4 3 1 3 1 5 3 4 1 5
## 19 1 1 1 1 1 1 1 1 1 1
## 20 7 7 7 7 7 7 7 7 7 7
## HARE BEAVER PORCUPIN MARMOT CSQUIRR REDSQUIR SM_DEN BB_DENSITY BLKBR_BUFB
## 1 4 1 1 1 1 1 4 3 1
## 2 4 1 4 3 5 3 3 3 4
## 3 3 1 4 1 3 3 4 3 1
## 4 4 1 1 1 3 4 4 4 4
## 5 3 1 1 4 5 3 3 3 1
## 6 5 1 3 1 4 4 5 3 5
## 7 4 1 3 1 4 3 4 4 4
## 8 5 1 4 1 5 4 3 4 5
## 9 3 1 4 1 3 3 4 3 1
## 10 4 5 5 1 4 3 4 5 4
## 11 5 4 5 1 4 4 5 4 6
## 12 4 1 1 1 1 1 4 3 1
## 13 7 7 7 7 7 7 7 7 1
## 14 3 1 4 1 3 3 4 3 1
## 15 4 1 4 1 5 4 3 5 1
## 16 1 1 1 1 1 1 5 3 3
## 17 4 4 1 1 3 5 5 5 4
## 18 3 1 1 4 5 3 3 3 1
## 19 1 1 1 1 1 1 1 3 1
## 20 7 7 7 7 7 7 7 7 1
## G_BEAR_A
## 1 3
## 2 3
## 3 3
## 4 4
## 5 3
## 6 4
## 7 5
## 8 4
## 9 3
## 10 3
## 11 4
## 12 3
## 13 7
## 14 3
## 15 3
## 16 0
## 17 4
## 18 3
## 19 0
## 20 7# look at the projection of the shapefile (note the use of "@" instead of "$")
elc_habitat@proj4string@projargs## [1] "+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs"wolfyht@proj4string@projargs## [1] "+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs"Sometimes we will need to change the projection of the shapefile to another projection (package rgdal). To change the projection we need to use the approrpiate epsg (Geodetic Parameter Dataset) code that represents your new projection.
Find a list of these spatial reference codes here: http://spatialreference.org/ref/epsg/
elc_habitat <- spTransform(elc_habitat, CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"))## Warning: PROJ support is provided by the sf and terra packages among others# check new projection in geographic coordinate system WGS84
elc_habitat@proj4string@projargs## [1] "+proj=longlat +datum=NAD83 +no_defs"# reset projection back to what is was previously
elc_habitat <- spTransform(elc_habitat, CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"))
# look at the projection of the shapefile
elc_habitat@proj4string@projargs## [1] "+proj=longlat +datum=NAD83 +no_defs"# another way to change it back to previous
elc_habitat <- spTransform(elc_habitat, CRS("+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
#write shapefile to files
writeOGR(elc_habitat,".","elc_habitat_NEW",driver="ESRI Shapefile", overwrite_layer=TRUE)## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others5.2 Extents
Spatial data objects have a spatial extent, in geographic space, that they occupy. To look at the spatial extent of the shapefile,we use ?extent()
extent(elc_habitat)## class : Extent
## xmin : 443680.6
## xmax : 627727.9
## ymin : 5618405
## ymax : 5789236extent(wolfyht)## class : Extent
## xmin : 555853
## xmax : 605389
## ymin : 5656997
## ymax : 5741316We will learn later how to harmonize datasets of different spatial extents.
6 Working with Rasters in R
par(mfrow= c(1,1)) ## reset graphical parameters
# reading in raster files (raster package)
deer_w<-raster("Data/deer_w2.tif")
moose_w<-raster("Data/moose_w2.tif") ## missing moose
elk_w<-raster("Data/elk_w2.tif")
sheep_w<-raster("Data/sheep_w2.tif") ## missing sheep
goat_w<-raster("Data/goat_w2.tif")
wolf_w<-raster("Data/wolf_w2.tif")#
elevation2<-raster("Data/Elevation2.tif") #resampled
disthumanaccess2<-raster("Data/DistFromHumanAccess2.tif") #resampled
# make a very basic plot of raster
plot(deer_w)
# look at the class of the raster
class(deer_w)## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"# look at basic raster summary
deer_w## class : RasterLayer
## dimensions : 5694, 6892, 39243048 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 443680.6, 650440.6, 5618416, 5789236 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=11 +datum=NAD83 +units=m +no_defs
## source : deer_w2.tif
## names : deer_w2
## values : 0, 7 (min, max)# look at raster data
deer_w@data## An object of class ".SingleLayerData"
## Slot "values":
## logical(0)
##
## Slot "offset":
## [1] 0
##
## Slot "gain":
## [1] 1
##
## Slot "inmemory":
## [1] FALSE
##
## Slot "fromdisk":
## [1] TRUE
##
## Slot "isfactor":
## [1] FALSE
##
## Slot "attributes":
## list()
##
## Slot "haveminmax":
## [1] TRUE
##
## Slot "min":
## [1] 0
##
## Slot "max":
## [1] 7
##
## Slot "band":
## [1] 1
##
## Slot "unit":
## [1] ""
##
## Slot "names":
## [1] "deer_w2"# look at structure of raster
str(deer_w)## Formal class 'RasterLayer' [package "raster"] with 12 slots
## ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots
## .. .. ..@ name : chr "C:\\Users\\Administrator.KJLWS11\\Documents\\Lab1\\Data\\deer_w2.tif"
## .. .. ..@ datanotation: chr "FLT8S"
## .. .. ..@ byteorder : chr "little"
## .. .. ..@ nodatavalue : num -Inf
## .. .. ..@ NAchanged : logi FALSE
## .. .. ..@ nbands : int 1
## .. .. ..@ bandorder : chr "BIL"
## .. .. ..@ offset : int 0
## .. .. ..@ toptobottom : logi TRUE
## .. .. ..@ blockrows : Named int 1
## .. .. .. ..- attr(*, "names")= chr "rows"
## .. .. ..@ blockcols : Named int 6892
## .. .. .. ..- attr(*, "names")= chr "cols"
## .. .. ..@ driver : chr "gdal"
## .. .. ..@ open : logi FALSE
## ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots
## .. .. ..@ values : logi(0)
## .. .. ..@ offset : num 0
## .. .. ..@ gain : num 1
## .. .. ..@ inmemory : logi FALSE
## .. .. ..@ fromdisk : logi TRUE
## .. .. ..@ isfactor : logi FALSE
## .. .. ..@ attributes: list()
## .. .. ..@ haveminmax: logi TRUE
## .. .. ..@ min : num 0
## .. .. ..@ max : num 7
## .. .. ..@ band : int 1
## .. .. ..@ unit : chr ""
## .. .. ..@ names : chr "deer_w2"
## ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots
## .. .. ..@ type : chr(0)
## .. .. ..@ values : logi(0)
## .. .. ..@ color : logi(0)
## .. .. ..@ names : logi(0)
## .. .. ..@ colortable: logi(0)
## ..@ title : chr(0)
## ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots
## .. .. ..@ xmin: num 443681
## .. .. ..@ xmax: num 650441
## .. .. ..@ ymin: num 5618416
## .. .. ..@ ymax: num 5789236
## ..@ rotated : logi FALSE
## ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
## .. .. ..@ geotrans: num(0)
## .. .. ..@ transfun:function ()
## ..@ ncols : int 6892
## ..@ nrows : int 5694
## ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs"
## .. .. ..$ comment: chr "PROJCRS[\"unknown\",\n BASEGEOGCRS[\"unknown\",\n DATUM[\"North American Datum 1983\",\n E"| __truncated__
## ..@ history : list()
## ..@ z : list()# look at the projection of the raster (note the use of "@" instead of "$")
deer_w@crs@projargs## [1] "+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs"# look at the spatial extent of the raster
extent(deer_w)## class : Extent
## xmin : 443680.6
## xmax : 650440.6
## ymin : 5618416
## ymax : 5789236# look at the resolution of the raster
res(deer_w)## [1] 30 30And if we need to change the projection of the raster to another projection (package rgdal) - note for the Tutorial I’ve masked these out # because they take a LONG time.
#deer_w2 <- projectRaster(deer_w, crs="+init=epsg:4326")
## or try this the same way we did above
#deer_w2 <- projectRaster(deer_w, crs="+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0")
#check projection of the raster
deer_w@crs@projargs## [1] "+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs"#change it back to what it was using another raster layers projection
#deer_w <- projectRaster(deer_w2, wolf_w)
#check projection of the raster
deer_w@crs@projargs## [1] "+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs"6.1 Creating a Raster Stack
One of the GREAT things about doing GIS in R is the ability to create a consistent raster stack of rasters of the same projection, extent, etc. LEts create a raster stack!
all.rasters<-stack(deer_w, moose_w, elk_w, sheep_w, goat_w, elevation2, disthumanaccess2)
plot(all.rasters)
#check class
class(all.rasters)## [1] "RasterStack"
## attr(,"package")
## [1] "raster"We can also sometimes need to export raster files to other formats such as a GeoTiff.
#learn more about
#?writeRaster()
writeRaster(deer_w, "Output/new_deer.tiff", "GTiff", overwrite=TRUE)7 Mapping with Mapview
Learn about mapview here: https://r-spatial.github.io/mapview/articles/articles/mapview_02-advanced.html
And also learn about different map types here http://leaflet-extras.github.io/leaflet-providers/preview/
Note, I had unexplained errors with the latest version of Mapview that were creating all sorts of problems both viewing the wolf telemetry points, and, rendering the .html through R markdown. I found the solution here : https://github.com/r-spatial/mapview/issues/312
install.packages("remotes")## package 'remotes' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\Administrator.KJLWS11\AppData\Local\Temp\Rtmp84rxPr\downloaded_packageslibrary(remotes)
remotes::install_github("r-spatial/mapview")## Downloading GitHub repo r-spatial/mapview@HEAD## cli (3.5.0 -> 3.6.0 ) [CRAN]
## raster (3.6-11 -> 3.6-13) [CRAN]
## htmlwidgets (1.6.0 -> 1.6.1 ) [CRAN]
## svglite (2.1.0 -> 2.1.1 ) [CRAN]## Installing 4 packages: cli, raster, htmlwidgets, svglite## Warning: package 'raster' is in use and will not be installed##
## There are binary versions available but the source versions are later:
## binary source needs_compilation
## cli 3.5.0 3.6.0 TRUE
## svglite 2.1.0 2.1.1 TRUE
##
## Binaries will be installed
## package 'cli' successfully unpacked and MD5 sums checked## Warning: cannot remove prior installation of package 'cli'## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying
## C:\Program Files\R\R-4.2.2\library\00LOCK\cli\libs\x64\cli.dll to C:\Program
## Files\R\R-4.2.2\library\cli\libs\x64\cli.dll: Permission denied## Warning: restored 'cli'## package 'htmlwidgets' successfully unpacked and MD5 sums checked
## package 'svglite' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\Administrator.KJLWS11\AppData\Local\Temp\Rtmp84rxPr\downloaded_packages## Running `R CMD build`...## * checking for file 'C:\Users\Administrator.KJLWS11\AppData\Local\Temp\Rtmp84rxPr\remotes7bc8424375d2\r-spatial-mapview-64f865a/DESCRIPTION' ... OK
## * preparing 'mapview':
## * checking DESCRIPTION meta-information ... OK
## * checking for LF line-endings in source and make files and shell scripts
## * checking for empty or unneeded directories
## * building 'mapview_2.11.0.9005.tar.gz'library(mapview)##
## Attaching package: 'mapview'## The following object is masked _by_ '.GlobalEnv':
##
## trailsmapView(wolfyht, layer.name = , zcol="NAME", native.crs = FALSE, legend = TRUE, cex=5, lwd=2, map.type = "Esri.DeLorme") + mcp2mapView(wolfyht, zcol="NAME", legend = TRUE, cex=5, lwd=2, map.type = "Esri.WorldImagery") + mcp2## Open Topo Map
mapView(wolfyht, zcol="NAME", legend = TRUE, cex=5, lwd=2, map.type = "OpenTopoMap") + mcp2## Esri World Topo Map
mapView(wolfyht, zcol="NAME", legend = TRUE, cex=5, lwd=2, map.type = "Esri.WorldTopoMap") + mcp2